function [pvSurplusesConsNum,pvRiskSharingConsNum,logValueNoBonds,sigmaTotalWealthNoBonds,aggrPriceOfRiskNoBonds] ...
  = solveModelAdditionalPdes(xGrid,xDrift,xVolatility,modelSolution,params,algParams)
%SOLVEMODELADDITIONALPDES Solves additional pdes not solved by solveModelPde()
%   Specifically, this function determines the bond value decomposition
%   into a surplus and a risk sharing (service flow) term and it determines
%   the value function in a model without safe assets. Input arguments are
%   analogous to solveModelPde(). The additional parameter modelSolution is
%   the output structure from that function.

  [theta,logValue,sigmaI,totalWealth,adjBondGrowth,aggrPriceOfRisk,qB,a,govSpending] = aux.unpack(modelSolution,'theta','logValue','sigmaI','totalWealth','adjBondGrowth','aggrPriceOfRisk','qB','a','govSpending');

  % unpack params
  [sigma,chi,gamma,rho,phi,iota0,delta,a0,g0,adjBondGrowth0,alphaA,alphaBondGrowth] = ...
    aux.unpack(params,'sigma','chi','gamma','rho','phi','iota0','delta','a0','g0','adjBondGrowth0','alphaA','alphaBondGrowth');
  % unpack algParams
  [nSolutionSteps,dt,convCrit] = aux.unpack(algParams,'nSolutionSteps','dt','convCrit');
  
  % prepare solution iterative loop
  differentiator = finiteDifferentiator(xGrid,3);


  diffusion = 1/2*xVolatility.^2;
  diffusion([1,end]) = 0;
  advection = xDrift;
  discountingPv = -rho*ones(size(xGrid));
  discountingLogValue = discountingPv;
  
  
  sigmaValue = differentiator.d1(logValue).*xVolatility;
  
  % initial guesses
  pvRiskSharing = theta;
  pvSurpluses = zeros(size(xGrid));
  logValueNoBonds = logValue;

  % main iterative loop to solve PDEs
  for i=1:nSolutionSteps
    % compute volatilities
    volatilityPvSurpluses = differentiator.d1(pvSurpluses).*xVolatility;
    volatilityPvRiskSharing = differentiator.d1(pvRiskSharing).*xVolatility;
    sigmaValueNoBonds = differentiator.d1(logValueNoBonds).*xVolatility;
    % compute stuff for logValueNoBonds
    iotaNoBonds = iota0 + ((a-govSpending-iota0)-rho)./(1+phi*rho);
    growthRateNoBonds = investmentTechnology(iotaNoBonds,phi,iota0) - delta;
    totalWealthNoBonds = (1 + phi*(a-govSpending-iota0))./(1+phi*rho);
    consCapRatioNoBonds = rho*totalWealthNoBonds;
    idiosyncraticRiskPremiumNoBonds = gamma*chi^2.*sigmaI.^2;
    inhomogeneityLogValueNoBonds = -(gamma-1)*(rho*log(consCapRatioNoBonds) + growthRateNoBonds - idiosyncraticRiskPremiumNoBonds/2) + sigmaValueNoBonds.^2/2;

    % update pv surpluses
    surpluses = - adjBondGrowth.*theta;
    inhomogeneitySurpluses = surpluses + sigmaValue.*volatilityPvSurpluses;
    pvSurplusesLong = linearParabolicPDE_timeStep([xGrid(1)-1;xGrid;xGrid(end)+1],[0;pvSurpluses;0],...
      dt,[NaN;diffusion;NaN],[NaN;advection;NaN],[NaN;discountingPv;NaN],[NaN;inhomogeneitySurpluses;NaN],0,0);
    pvSurplusesNew = pvSurplusesLong(2:end-1);
    % update pv risk sharing
    riskSharingBenefit = gamma*(1-theta).^2.*chi^2.*sigmaI.^2.*theta;
    inhomogeneityRiskSharing = riskSharingBenefit + sigmaValue.*volatilityPvRiskSharing;
    pvRiskSharingLong = linearParabolicPDE_timeStep([xGrid(1)-1;xGrid;xGrid(end)+1],[0;pvRiskSharing;0],...
      dt,[NaN;diffusion;NaN],[NaN;advection;NaN],[NaN;discountingPv;NaN],[NaN;inhomogeneityRiskSharing;NaN],0,0);
    pvRiskSharingNew = pvRiskSharingLong(2:end-1);
    
    % update logValueNoBonds
    logValueNoBondsLong = linearParabolicPDE_timeStep([xGrid(1)-1;xGrid;xGrid(end)+1],[0;logValueNoBonds;0],...
      dt,[NaN;diffusion;NaN],[NaN;advection;NaN],[NaN;discountingLogValue;NaN],[NaN;inhomogeneityLogValueNoBonds;NaN],0,0);
    logValueNoBondsNew = logValueNoBondsLong(2:end-1);

    % compute convergence metric
    relChangeSurpluses = (pvSurplusesNew-pvSurpluses)./(abs(pvSurplusesNew)+abs(pvSurpluses))*2;
    relChangeRiskSharing = (pvRiskSharingNew-pvRiskSharing)./(abs(pvRiskSharingNew)+abs(pvRiskSharing))*2;
    relChangeLogValueNoBonds = (logValueNoBondsNew-logValueNoBonds)./(abs(logValueNoBondsNew)+abs(logValueNoBonds))*2;
    maxRelChangeSurpluses = max(abs(relChangeSurpluses));
    maxRelChangeRiskSharing = max(abs(relChangeRiskSharing));
    maxRelChangeLogValueNoBonds = max(abs(relChangeLogValueNoBonds));
    maxRelChange = max([maxRelChangeSurpluses,maxRelChangeRiskSharing,maxRelChangeLogValueNoBonds]);

    % transfer variables
    pvSurpluses = pvSurplusesNew;
    pvRiskSharing = pvRiskSharingNew;
    logValueNoBonds = logValueNoBondsNew;

    if maxRelChange < convCrit
      break;
    end
  end
  
  pvRiskSharingConsNum = pvRiskSharing.*totalWealth;
  pvSurplusesConsNum = pvSurpluses.*totalWealth;
  
  % model without bonds
  totalWealthNoBonds = (1 + phi*(a-govSpending-iota0))./(1 + phi*rho);
  sigmaTotalWealthNoBonds = differentiator.d1(log(totalWealthNoBonds)).*xVolatility;
  sigmaValueNoBonds = differentiator.d1(logValueNoBonds).*xVolatility;
  aggrPriceOfRiskNoBonds = -sigmaValueNoBonds + sigmaTotalWealthNoBonds;
  
end

